Laboratory Exercise 10
Multilevel modeling 2: Random intercept and random slopes models
Part 1 - Data importing and tidying
Starting up
As per last week, download and unzip today’s data into a new folder and open the file “tutorial_10.Rproj”. This will open up R and R studio and make sure that you don’t have to worry about file-paths.
Setting up R to run our analyses
Next, try running the following code by copying and pasting (Ctrl+C, Ctrl+V) it into your R session’s console.
# The following code functions ensure that the packages (like sub-programs that R uses to do things) are installed and attached
for(package in c('lme4', 'lmerTest', 'tidyverse', 'foreign', 'readr', 'sjstats')) {
if(!require(package, character.only=T, quietly=T)) {
install.packages(package)
library(package, character.only=T)
}
}The data
The data for this week’s exercise comes is the same as last weeks - it is an organizational psychology data-set described in Klein et al (2000), and is stored as siop_merged.csv.
Variables
The file siop_merge.csv contains 750 employee-level observations nested within 50 work-groups. Ignoring the variable group ID (grpid), there are eight standardized variables, where a higher score indicates a higher level (e.g., higher pay or more negative leadership behaviors). These data were collected by individual survey of the employees, so measure individual level perceptions, except for the variable “physen”, which is the group level “physical work environment”.
Table 1. The variable names and variables included in siop_merge.csv.
The first thing we need to do is to read in the data. Copy and paste the following into your R console to load today’s data.
siop_merged <- read.csv("siop_merged.csv")
The data
The first variables we’ll be dealing with today are “Job satisfaction”, “Positive affect” and “physical work environment (group level)”. As always its a good idea to plot your data. This data has three variables, so you can either plot it with facets, or use a 3D plot. Click and drag this one to have a look at your data;
Figure 1. A 3D scatterplot of Job satisfaction, positive affect and group Physical environment, color indicates group membership.
Question 1.
What is the relationship between positive affect and job satisfaction? Is there a clear directional relationship?
Part 2 - Random intercepts model with a predictor at the group level
Let’s fit a random intercept model with one individual level predictor - we’ll start with positive affect - and one group level predictor - physical work environment.
This model with random intercepts for group and a single level 1 predictor (positive affect) may be written:
\(JobSat_{ij} = \gamma_{00} + \gamma_{10} PosAff_{ij}+ \gamma_{01}PhysEn_{0j} + u_{0j} + \epsilon_{ij}\)
This means each person’s job satisfaction is equal to the overall intercept ( \(\gamma_{00}\) ), a regression term times their positive affect ( \(\gamma_{10}(PosAff_{ij})\) ), a regression term times the group level predictor ( \(\gamma_{01}PhysEn_{0j}\) ), a random effect for their group (\(u_{0j}\)), plus an error term (\(\epsilon_{ij}\)) (i.e., the residual, as the MLM will not perfectly predict their score). You can run this model with the following code:
PA_PE_model <- lmer("jobsat ~ 1 + posaff + physen + (1 | grpid)", data = siop_merged)Which should be read as:
| Terms | Meaning |
|---|---|
| jobsat ~ | job satisfaction is predicted by |
| 1 + | An intercept parameter which is the same for everyone |
| posaff | a regression coefficient times each person’s positive affect |
| physen | a regression coefficient times each group’s physical enviroment |
| (1 | grpid) | A random intercept for each group |
Remember summary(PA_PE_model) is useful in getting a summary of your output out!
summary(PA_PE_model)This model says that the relationship between positive affect and job satisfaction is the same in every group, but that the baseline - the intercept - can vary by group, and that we can predict the group intercepts using the group’s score on physical environment. To visualize this, we can plot job satisfaction by positive affect and draw a regression line for each group. The model does not actually specify an intercept for each group, but rather assumes that the intercepts come from a normal distribution with a mean of the overall intercept and a variance as specified in the random effects output above. The output we will often be interested is the variance of these groups’ intercepts as opposed to their individual values.
Figure 2
Figure 2. Job satisfaction by positive affect, lines show the estimated relationship by group, color indicates group membership.
Question 2.
What is the relationship between positive affect and job satisfaction? Is there a clear directional relationship?
Question 3.
Looking at the following equation, can you identify what model ouput corresponds to \(u_{0j}\)?
\(JobSat_{ij} = \gamma_{00} + \gamma_{10} PosAff_{ij}+ \gamma_{01}PhysEn_{0j} + u_{0j} + \epsilon_{ij}\)
Hint: This may be a trick question.
Question 4.
How much variance is accounted for by the grouping factor (grpid)?
Hint: Calculate the ICC for this model.
Part 3 - Random intercepts and random slopes
Now we will include a random slope for the association between positive affect and job satisfaction as well as a group level predictor for the intercept, but without other predictors for the slope. This means that positive affect becomes both a fixed and a random effect.
This model with random intercepts and a random slope for group, a single group level predictor (physical work environment), and a single level 1 predictor (positive affect) may be written:
\(JobSat_{ij} = \gamma_{00} + \gamma_{10} PosAff_{ij} + \gamma_{01} PhysEn_{0j} + u_{1j}PosAff_{ij} + u_{0j} + \epsilon_{ij}\)
This means each person’s job satisfaction is equal to the overall intercept ( \(\gamma_{00}\) ), a regression term times their positive affect \(\gamma_{10} (PosAff_{ij})\), a regression term times the group level predictor \(\gamma_{01}PhysEn_{0j}\), a random effect for positive affect ( \(u_{1j}PosAff_{ij}\) ), a random effect for their group (\(u_{0j}\)), plus an error term (\(\epsilon_{ij}\)) (i.e., the residual, as the MLM will not perfectly predict their score).
RS_model <- lmer("jobsat ~ 1 + posaff + physen + (1 + posaff | grpid)", data = siop_merged)Which should be read as:
| Terms | Meaning |
|---|---|
| jobsat ~ | job satisfaction is predicted by |
| 1 + | An intercept parameter which is the same for everyone |
| posaff | a regression coefficient times each person’s positive affect |
| physen | a regression coefficient times each group’s physical enviroment |
| (1 + posaff | grpid) | A random intercept and a random slope for the relationship between positive affect and job satisfaction for each group |
Remember summary(RS_model) is useful in getting a summary of your output out!
summary(RS_model)That this model says that the relationship between positive affect and job satisfaction changes in different groups by some estimated amount, and that we can predict the groups’ true mean scores using their group’s physical environment. To visualize this, we can plot job satisfaction by positive affect and draw a regression line for each group. The regression line slopes and intercepts change by group, because we now have random slopes and intercepts, allowing the relationship between positive affect and job satisfaction to change by group. Conceptually, this model is better thought of as assuming that the group slopes and intercepts come from a normal distribution some overall mean (i.e., the fixed effects for the intercepts and slopes) and an estimated variance. The output we’re often interested in above and beyond the fixed effects is the variance of these groups’ intercepts, and the variance of the groups slopes.
Figure 3
Figure 3. Job satisfaction by positive affect, lines show the estimated relationship by group, color indicates group.
Q
What does it mean to include random slopes in the model? What are we saying about the relationship between job satisfaction and positive affect? Is that claim plausible or even likey?
Part 4 - Random intercepts and random slopes with a group level predictor
Let’s fit a random intercept and slopes model with one individual level predictor - again positive affect - and one group level predictor - physical work environment, and one group level predictor of the random effects.
This model with random slopes and intercepts for group, a single level 1 predictor (positive affect), and an interaction between workload and physical environment may be written:
\(JobSat_{ij} = \gamma_{00} + \gamma_{10} PosAff_{ij}+ \gamma_{01}PhysEn_{0J} + u_{0j} + \gamma_{11}(PosAff_{ij}\times PhysEn_{0J}) + \epsilon_{ij}\)
This means each person’s job satisfaction is equal to the overall intercept ( \(\gamma_{00}\) ), a regression term times their positive affect \(\gamma_{10} (PosAff_{ij})\), a regression term times the group level predictor \(\gamma_{01}PhysEn + u_{0j}\), a random effect for their group (\(u_{0j}\)), plus an error term for the individual (\(\epsilon_{ij}\)).
RS_model_2 <- lmer("jobsat ~ 1 + physen + posaff + posaff * physen + (1 + posaff | grpid)", data = siop_merged)Which should be read as:
| Terms | Meaning |
|---|---|
| jobsat ~ | job satisfaction is predicted by |
| 1 + | An intercept parameter which is the same for everyone |
| posaff | a regression coefficient times each person’s positive affect |
| physen | a regression coefficient times each group’s physical enviroment |
| physen:posaff | a regression coefficient times the product of physical enviroment times positive affect |
| (1 + posaff | grpid) | A random intercept and a random slope for the relationship between positive affect and job satisfaction for each group |
Remember summary(PA_PE_model) is useful in getting a summary of your output out!
summary(RS_model_2)That this model says that the relationship between positive affect and job satisfaction changes in different groups by some estimated amount, and that we can predict the groups’ true mean scores using their group’s physical environment, and that the relationship between positive affect and job satisfaction is influenced by the physical environment.
To visualize this, we can again plot job satisfaction by positive affect and draw a regression line for each group. The regression line slopes and intercepts change by group, because we now have random slopes and intercepts, allowing the relationship between positive affect and job satisfaction to change by group.
Figure 4
Figure 4. Job satisfaction by positive affect, lines show the estimated relationship by group, color indicates group.